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Global linear stability analysis of a self-similar solution describing the interaction of a relativistic shell with an ambient medium is per- 
formed. The solution is shown to be unstable to convective Rayleigh- Taylor modes having angular scales smaller than the causality scale. 
Longer wavelength modes are stable and decay with time. For modes of sufficiently large spherical harmonic degree I the dimensionless 
growth rate scales as yJl/V, where F is the Lorentz factor of the shell. The instability commences at the contact interface separating 
the shocked ejecta and shocked ambient gas and propagates to the shocks. The reverse shock front responds promptly to the instability 
and exhibits rapidly growing distortions at early times. Propagation to the forward shock is slower, and it is anticipated that the region 
near the contact will become fully turbulent before the instability is communicated to the forward shock. The non-universality of the 
Blandford-McKee blast wave solution suggests that turbulence generated by the instability in the shocked ambient medium may decay 
slowly with time and may be the origin of magnetic fields over a long portion of the blast wave evolution. It is also speculated that the 
instability may affect the emission from the shocked ejecta in the early post-prompt phase of GRBs. 
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1 Introduction 

Relativistic shock waves is a common phenomenon in astrophysics. They form when a relativistic flow 
ejected by a compact central engine interacts with the surrounding medium or, in case of an intermittent 
source, as a result of steepening of waves produced in the outflow itself. The broadband emission observed in 
blazars, micro-quasars and gamma-ray bursts (GRBs) is produced behind those shocks and is an important 
diagnostic of the dissipation process. 

Despite an impressive progress in our understanding of relativistic shocks some outstanding problems 
remain open. Of particular interest is the relativistic blast wave that form in GRB explosions when the 
relativistic ejecta expelled by the source impact the circumburst medium. The afterglow emission observed 
in most long GRBs is most likely produced in the thin layer enclosed between the forward shock and 
the ejecta, and is an important diagnostic of the blast wave evolution and the conditions in the shocked 
layer. Although the simple blast wave model has been quite successful in explaining the late afterglow 
evolution, recent observational efforts revealed some features that require extension of the simple model: 
(i) Observations of the late afterglow emission indicate strong amplification of magnetic fields in the post 
shock region - by several orders of magnitudes larger than what can be achieved by compression of the 
ambient magnetic field. Despite recent efforts to investigate potential mechanisms by which magnetic 
fields can be generated or amplified in the vicinity of the shock, this issue remains unresolved, (ii) SWIFT 
observations during the early afterglow phase reveal strong deviation of the lightcurve at early times 
from that predicted by the simple blast wave model. Several ad hoc explanations have been offered, 
including prolonged activity of the central engine and evolution of microphysical parameters. However, the 
feasibility of these scenarios depends on poorly understood physics, and it remains to be demonstrated 
that they can be derived from first principles, (iii) In the fireball scenario commonly adopted, the naive 
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expectation has been that the crossing of the reverse shock should produce an observable optical flash. 
Despite considerable observational efforts, suc h flashes seem to be very rare. One plausible explanation 
is that the electa is magnetically dominated (|Levinson and Eichler 19931 . Lyutikov and BlandfordI 2003 , 
Giannios and Spruit 20051 . Giannios et al. 20081 ). Poynting flux dominated outflows have the advantage 



that they can naturally account for the ultra-relativistic Lorentz factors inferred. On the other hand, they 
are challenged by the rapid dissipation of magnetic energy that seems to be required. Furthermore, even 
if the flow is magnetically dominated some accumulation of baryon rich matter at the 'piston's' head is 
anticipated during the shock breakout phase, that may mimic effects of a hydrodynamic ejecta. 

In this paper we explore the stability of the double-shock system. Hydrodynamic instabilities can give 
rise to strong distortions of the structure that may generate turbulence, amplify magnetic fields, and affect 
the emission processes in the post-prompt phase. Such effects have been studied in the non-relativistic 
case in connection with young supernovae remnants (SNRs). In fact, the idea that the Rayleigh- Taylor 
(R-T) instability should play an important role in the deceleration of a non-relativistic ejecta dates back to 
Gull (I1973I). who perf ormed ID simulations of young SNRs that incorporate a simple model of convection. 
IChevalier et al. ( 19921 ) later performed a global linear stability analysis of a self-similar solution describing 
the interaction of non-relativistic ejecta with an ambient medium and found that it is subject to a convective 
instability. They analyzed self-similar perturbations and showed that the flow is unstable for modes having 
angular scales smaller than some critical value. The convective growth rate was found to be largest at the 
contact discontinuity surface and to increase with increasing / number of the eigenmodes. They also 
performed 2D hydro dynamical simulations that verified the linear results and enabled them to study the 
nonlinear evolution of the instability. The simulation exhibits rapid growth of fingers from the contact 
interface that saturates, in the nonlinear state, by the Kelvin-Helmholtz instability. Str ong distortions of 
the co ntact and the reverse shock was observed with little effect on the forward shock. lJun and Norman 
(jl996l ) performed 2 and 3D MHD simulations of the instability to study the evolution of magnetic fields in 



the c onvection zone. They confirmed the rapid growth of small scale structure reported in (jChevalier et al. 
1993), and in addition found strong amplification of ambient magnetic fields in the turbulent flow around 



R-T fingers. On average, the magnetic field energy density reaches about 0.5% of the energy density of the 
turbulence, but it could well be that t he magnetic field ampl ificat ion was limited by nurn erical resolution 
in their simulations. The simulations of Chevalier et al. (|l992l ^ and lJun and NormanI (Il99fil ^ support earlier 
ideas, that the clumpy shell structure observed in young (pre-Sedov stage) SNRs such as Tycho, Kepler 
and Gas A is due to the R-T and K-H instability. 



In this paper we extend the linear stability analysis of lGhevalier et al\ (1992 ) into the relativistic regime. 



A preliminary account of the model and results is presented in (iLevinson 



20091 ). We find that denser ejecta 



sweeping a lighter ambient gas are subject to the R-T instability also in the relativistic case. The reason is 
that in the rest frame of the decelerating contact there is an effective gravitational force which is directed 
outwards, and so in this frame the denser ejecta is 'on top' o f the lighter ambient gas. The stability of 
a double-shock system has been investigated by Wang et al. ( 20021 ) using the thin shell approximation. 
However, this study is limited to large scale mod es and neglects pressure gradients and, therefore, pre- 
cludes the convective instability. Thompson ( 20061 ) pointed out that a magnetized photon-rich shell that 
propagates through a dense Wolf-Rayet wind may be subject to the R-T instability. Using heuristic ar- 
guments he examined the conditions under which the instability develops and estimat ed th e growth rate. 
The scaling of the growth rate found below is consistent with his result. Gruzinov (l200d) performed a 
linear stability analysis of a Blandford-McKee (BMK) blast wave solution (jBlandford and McKeel[l976l 'l. 
and found that the BMK solution is stable but non-universal, in the sense that some modes decay very 
slowly as the system evolves. Furthermore, the onset of oscillations of an eigenmode of order I has been 
seen in the simulation once the Lorentz factor evolved to F < /. The conclusion drawn based on Gruzinov's 
findings is that distortion of the shock front at early times may cause significant oscillations during a 
large portion of its evolution. If the amplitude of these oscillations is sufficiently large, and if the same 
behavi or holds in the nonlinear regime then this can lead to generat ion of vorticity in the post shock 
region ( Goodman and MacFadven 20081 . Milosavlievic and Nakai] 12007 1. and the consequent amplification 
of magnetic fields, as demonstrated recently by Zhang et al. ( 2009l V 
The plan of the paper is as follows: In section [2] we derive the basic equations in a general form. In 
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section [3] a class of self-similar solutions for the double-shock structure, obtained originally by Nakamura 
and Shigeyama (2006), is reviewed. These are employed as the unperturbed solutions for our analysis. The 
linear perturbation analysis of these solutions is presented in section HI The implications for gamma-ray 
bursts are discussed in section [H We conclude in section [6l Detailed derivation of main results is given in 
the appendices. 



2 Basic equations 

Consider an unmagnetized fluid, and let p, p, h and denote its proper density, pressure, dimensionless 
specific enthalpy, and 4-velocity, respectively. The stress-energy tensor then takes the form 

T^"" = phu^'u" - g^"'p, (1) 

where g'^'^ is the metric tensor. Neglecting radiative losses, the dynamics of the flow is governed by mass 
and energy-momentum conservation: 

d^{pu>') = 0, S^T^'^ = 0. (2a,b) 

Using ([2^) the different components of ([2)3) reduce to 

dln7 dp 

d 
dt 

Pl-Jlihl^r) + ^TP = 0. (3c) 
at 

Here 7 = is the Lorentz factor of the fluid, vt is the tangential component of the 3-velocity, which 
we express as v = f r r + vj-, 7 is the adiabatic index, W = ph^'^, d/dt = (ii^/ii°)9^ is the convective 
derivative, and 

+ (4) 
r od r sm a ocp 

If the flow passes through a discontinuous shock front, then the solutions of the flow equations in the 
upstream and downstream regions are to be matched at the shock surface, which is deflned by the equation 
ip{x'^) = r — R{t, 0, (p) = 0. Integration of (l2^,b) across the surface lead to the jump conditions 

[pu^]n^ = 0, [T'^n riu = 0, (5a,b) 

where the square brackets denote the difference of the enclosed quantity across the shock front, and 



^ln(p/p^)=0, (3b) 



(6) 



is a 4- vector normal to the shock front. 
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Figure 1. Schematic representation of the double-shock system. There are three characteristic surfaces: a forward shock propagating in 
the ambient medium, a reverse shock sweeping the ejecta, and a contact discontinuity separating the shocked ejecta and the shocked 
ambient medium. The Lorentz factors of the three surfaces, measured with respect to the unshocked ambient medium, are indicated. 
Quantities in the shocked ambient medium (region 1) and shocked ejecta (region 2) are denoted by subscripts 1 and 2, respectively. 



3 Unperturbed solutions 

The u nperturbed solution invoked here is the self-similar solution derived by Nakamura and Shigeyama 
(|2006l l. It is reviewed here to set up the notation and to introduce some aspects that are important for 



the stability analysis. Nakamura and Shigeyama considered a freely expanding ejecta interacting with an 
ambient medium having a density profile pi = br~^. The freely expanding ejecta is characterized by a 
velocity Ve = r /t aX time t after the explosion, and a proper density profile 

(7) 



where 7e = — v1 is the corresponding Lorentz factor (it can be readily seen that the continuity 

equation ([2^) is satisfied for this choice of pe and Ve). 

The system under consideration is shown schematically in figure [TJ The subscript 1 refers to the shocked 
ambient medium and 2 to the shocked ejecta. The Lorentz factors of the forward shock, reverse shock and 
the contact discontinuity are denoted by ri(t), T2{t) and Tc{t)^ respectively. Self-similarity requires that 
they all have a similar time evolution, viz., = At~^, Tf = Bt~^, = Ct~"^, where A,B,C and m 
are constants determined upon matchin g the solutions in regions 1 a nd 2 at the contact discontinuity. The 
similarity parameter can be defined as (jBlandford and McKeelll97d l 



X=[l+2im + l)Tl]{l-r/t). (8) 



The shocks and the contact are surfaces of constant X) since the velocity of a constant x surface is 
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given by 



^ = 1- 

dt 2ri 



1-:^. (9) 



Mt)=l (1-^ dt' = t---^, (10) 



we readily obtain Xi = 1; Xc = i^i/^c)^ = B /C > 1 and X2 = (ri/r2)^ = B/A> Xc- The trajectory of 
the reverse shock is 

— dt' = t — — 

V '^^J 2{m + L,L^ 

from which we obtain for the velocity of the ejecta crossing the shock: Ve{r2) = r2/t = 1 — l/[2(m + l)r2]. 
The corresponding Lorentz factor is thus given, to order 0(r2 by 

^l = {m + l)Tl (11a) 

and the density by 

~ t^^^ ~ A^/-^{m + l)-/2' 2 • Uit^j 

3.1 Shocked ambient medium 

We consider cases where the forward shock is ultra-relativistic. The specific enthalpy of the shocked ambient 
gas (region 1) is then approximated by hi = 4pi/pi. The jump conditions at the forward shock follow from 
([5^1, b) and ([6]) using R{t,9,(j)) = ri{t) = t[l — l/2(m + l)rf]. The self-similar variables for the Lorentz 
factor, pressure and density are defined as 

7? = l^lgix), (12a) 
p[ = P171 = 2p,Tlh{x), (12b) 
Pi = lp^r^fix)■ (12c) 

Under this choice the shock jump conditions imply ^(1) = /(I) = h{l) = 1. The equations obeyed by 
the variables g, f, h are outlined in Appendix lA.il Using ()12h ) the Lorentz factor of the contact surface 
can be written as Tc = Tiy/gc/2, where gc = g{Xc), which combined with the relation Xc = (ri/Tc)^ 
derived above yields Xc9c = 2. The solution is obtained by numerically integrating equations ()A.2b -c) from 
the forward shock front x = 1 to the contact discontinuity where gcXc = 2. The quantities gc, Xc are 
eigenvalues of the solution. 

3.2 Shocked ejecta 

The reverse shock cannot be considered ultra-relativistic in general and, therefore, a complete treatment 
is required. The specific enthalpy of the shocked ejecta is taken to be /12 = 1 -|- 7P2/[P2(7 — l)]i and we 
remind that 7 denotes the adiabatic index. Assuming the unshocked ejecta to be cold, the jump conditions 
at the reverse shock, obtained from (l5^,b) and dGj) using R{t,6,(j)) = r2(t), read 

PeleiVe " V2) = P2l2{v2 - V2), (13a) 
Pell {Ve -V2) = W2{V2-V2)+ P2 V2 , (13b) 
PellveiVe - V2) = W2V2{v2 " U2) + P2, (13c) 
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Figure 2. Profiles of the pressure, proper density and Lorentz factor of the unperturbed flow, for two different choices of parameters. 
The forward shock is located at x = 1- The location of the reverse shock (X2) and the contact discontinuity (xc) are indicated. 

here V2 = dr2/dt is the shock 3- velocity. Equations (I13b -c) can be solved by employing (jllh ) and recalling 
that W2 = P2^272- finds 



72 (X2) = q^l, 



P2iX2 



mqpeje 



{m + l){q - 1) 



mpe 



1 - Vq/im + 1) 



a{q - 1) + 2 

where q is the only positive solution of the equation 

jx^ + (2 — ^)Vm + — (2 — 7)2; — 7\/m + \x = 0. 



(14a) 
(14b) 

(14c) 



(15) 



Following Nakamura and Shigeyama we find it convenient to transform in region 2 to a new 

similarity parameter, 



<y = X/X2 = {1 + 2(m + l)r2}(l - r/t). 
The reverse shock is then located at o" = 1 and the contact at 

o-c = Xc/x2 = ri/r^. 

The self similar variables of the shocked ejecta, G, F, are then defined as 

7I = qTlG{a), 



P'2 



P2 



mqpe-fe 



(m+ l)(g- 1) 



Hia), 



mpe 



1 - ^/q/{m + 1) F{a), 



(16) 

(17) 

(18a) 
(18b) 

(18c) 



a(g - 1) + 2 L 

and satisfy F{1) = G{1) = H{1) = 1. From equations (fTTj) and <![L8h ) we obtain GcCTc = 1/q at the 
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contact discontinuity. The equations obeyed by these self-similar variables are derived in Appendix IA.2[ 
The solution in this region is obtained upon integration of (|A.4b -c) from the reverse shock cr = 1 to the 
contact GcCJc = 1/9- 



3.3 Conditions at the contact surface 

Two conditions at the contact discontinuity fix the constants A, 5, and m. One condition is that there be 
no flow across the contact interface. This implies Fc = 7ic = 720 from which we obtain 



^2=- = 29-, (19) 



where equations (|T2a ) and (|T8a ) have been employed. The second condition is pressure balance, viz., 
Pi(^)Xc) = P2{t, ac)- This condition yields two relations. The first one, 

6 -2k 

m = —, (20) 

comes from the requirement that pi and p2 have the same time dependence. The second one is implied by 
equations (fT^ l and (fTBb ): 



4g[7(g - 1) + 2(7 - 1)] GJc fni + 1 



5^1+n/2 3^(^ _ 1) g^p^ \ CTc J V \ m+1 



n/2 



l-J—^\ . (21) 



The solution described above is valid for —l<m<3 — k ( Blandford and McKed Imdi ). The case 
m = 3 — k corresponds to an adiabatic impulsive b last wave (for which = 1 so that the contact 
surface is undefined). From ()A.2b -cl it can be shown ( Blandford and McKed 19761 ) that near the contact 



discontinuity the density of the ambient medium behaves as /i oc (2— ^x)" S^ith^i = {m—k)/{m+3k— 12). 
Thus, within the range of parameters for which the solution is valid h diverges at the contact for m < k 
and vanishes for m > k. Likewise, from (lA.4b -c) we find that the density of the shocked ejecta behaves as 
H oc {qGa-l)-^^, where 6*2 = (1 - 7)(mn - 6)/(107 + 2mn - 12) = (7 - 1)(6 + n/c)/[57(n + 2) - 2nA; - 12]. 
Within the allowed range of parameters 62 is always positive, so that H always diverges. This is seen in 
figure [21 where solutions obtained for an ejecta interacting with a uniform density medium (right panel) 
and a stellar wind (left panel) are exhibited. 



4 Global linear stability analysis 

We now consider linear per turbations of the self- similar solution described above. As shown below, unlike 



in the non-relativistic case (jChevalier et a/.lll992l ) the relativistic perturbation equations do not admit self- 
similar solutions, with the exception of the spherical mode. The reason is the inherent coupling, via the 
Lorentz factor, of the radial and tangential velocity perturbations 5vr and 6vt- Thus, numerical integration 
of the time dependent equations is required. 

4.1 Perturbation equations 

The perturbed variables are taken to be 



P' = P'oit^r) + 6p{t,r,9,(j)), 

p = Po{t,r) +6p{t,r,9,(j)), 

V = vo{r)r + 5vr{r, 9, 0)r -|- 6vT{r, 9, 4>), 



(22a) 
(22b) 
(22c) 
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here the subscript zero denotes the unperturbed flow. To simphfy the notation we define = dt + fo<9r 
to be the convective derivative with respect to the unperturbed flow. The Lagrangian change of some 
quantity Q is then given to flrst order by dQ/dt = d^((5o + ^Q) + SvrdrQo- The linearized equations for 
the perturbations, obtained upon substitution of the expressions (|22b -c) into the hydrodynamic equations 
dS^-c) are 

dti^p'/p'o) + drSvr + (2/r + dr In p'Q)5vr + V(5vt) = 0, (23a) 



'^t\—-l^+l—}+ Svrdr HpoPoI = 0, (23b) 
\Po Po 70/ 

p[,dO(/io7o5vT) + Vt5p = 0, (23c) 

Wo [d'li^iSvr) + 5vrdr In 7o] + ,5WdO(In7o) 

+jid^t6p + ji6vr [27odtVo + drPo] - dt5p = 0. (23d) 

Next, we apply the perturbation equations (i23k -d) to regions 1 and 2 (see figured]), using the self-similar 
solutions in each region for the unperturbed quantities. 

4.1.1 Perturbed flow in region 1. We expand the perturbations in spherical harmonics and use x ^iid 
r = Int as new independent variables in place of r, t. The perturbations in region 1 are then expressed as 

5p'i = p[oip{r,x)YU0,(t^), (24a) 

5pi=Pwip{T,x)yifh{e,(t>): (24b) 
1 

JviT = Ct(t, x)rVTYirn{e, </>), (24d) 



5vir = ^iR{r, x)yirh{0, 0), (24c) 



with p'^Q andpio given by (|12b .c). respectively. When these expressions are substituted in equations (P5b -d) 
a set of first order hyperbolic PDEs for the dimensionless amplitudes is obtained: 

dric. = ^^{Ao^pd^ip. + B^pip}, (25) 

where the indices q, (3 run over R, P, p, T. The details are given in Appendix IB.li The coefficients Aa/3 and 
Ba/3, given explicitly in (|B.4b -cl and ()B.5h -g). are functions of the self-similarity coordinate X) but are 
independent of r. 



4.1.2 Perturbed flow in region 2. Likewise, in region 2 we transform to the coordinates a, r and define 



^f>2 = P2oVp{T,(^)Yifh{0,4>), (26a) 

Sp2 = P20f?p(T, cr)YirhiO, 4>), (26b) 

5v2r = ^^r]R{T,(T)Yirh{B,(t)), (26c) 

5v2T = i]T{T,a)rVTYifh{0,4>). (26d) 

The unperturbed parameters of the shocked ejecta, /?2q, P20 are given by equations (fTSb .c). The derivation 
of the perturbation equations for the dimensionless perturbations rja is presented in Appendix IB. 21 The 
resultant set of equations is 
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drr]a = T,f3{Caf3dxr]f3 + Dapijp}, (27) 
with the coefficients Cap{(y) and Da/s^cr) given exphcitly in (fR8k -d) and (lR9k-h). 

4.2 Boundary conditions 



Equations (|25|) . (|27j) are solved subject to boundary conditions imposed at the shock fronts and at the 
contact discontinuity. We ahow perturbations of the shock fronts and the contact surface of the form 

Sra{t,9,<P) = *-^Yi^{e,4>), (28) 

here a = 1, 2, c refers to the forward shock, reverse shock, and the contact discontinuity, respectively. The 
corresponding perturbation of the 3-velocity at these surfaces is 

5Va = r-^[dr5a + (m + l)5a]YirniO, 4>). (29) 

Now, the Lagrange perturbation of some fluid quantity Q at the perturbed position of surface a is given 
by AaQ = idrQo)Sra + 5Q; e.g., Aap' = {drPQ)6ra + 5p', etc. Defining v^^ = /vP and denoting by 
Uj^ = vP-^ + (5nj^ the perturbed normal of the forward (j = l)/reverse (j = 2) shock, equation ([5^) gives 
to first order 

[(p[) + A,pO« + ^jv^in^ + Km)] = 0, (30a) 

and dSb) 

[{Wo + A,WW^ + Ajv^iv^ + Ajv'') - (po + Ajp)g'"']{n% + dnj^) = 0. (30b) 

The forward shock is described by the equation ijji{x^) = r — ri{t) — 5ri{t,6,4>), from which we obtain, 
using ([6]), 

< = (-ril^i,ri,0), (31a) 



5ni^ = {-rl6Vi,rlVidVi,-riVT6n). (31b) 



The density, pressure, and Lorentz factor of the unshocked flow just upstream of the forward shock are 
Pi = 0, Pi = br~'^, and = 1. At the perturbed shock front Aipi = —kpi6ri/ri and can be neglected to 
the order to which we are working. Applying (j24b -d) to the forward shock and using ()30h ). (|31h .b) we get 

= d^6i + 2CR + {3m-3 + 2k)6i. (32a) 
Likewise, the transverse component of (pOb ) yields 

Ct = -di/rl (32b) 

and the other two components 

= 2drSi + (14 - 8m - 6A;) ^ , (32c) 

S,R = 2dr5i-2{m-3 + k)5i. (32d) 
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Finally, we eliminate 5i from equations (I32b .d) to get three boundary conditions at the forward shock 
(X=l): 

dr^P - = (m - 3 + k)Cp + ^{7 -Am- 3k)^R, (33a) 

^t = ^(^^Kp-C.], (33b) 
= -driTj^T) + 2Cr - (3m - 3 + 2k)TlCT. (33c) 
As a check note that for the imp ulsive BMK solution with m = 3 and k = equations (I33b .b) reduce to 



those derived in ( lGruzinovll200 0l) 



Id). 



The derivation of the boundary conditions at the reverse shock is far more involved. The details are 
presented in Appendix O One finds 

- fspdsp r, 1 - fsSdsp „^ , . 

dTr]R = --j f ^^ rVP--^ r . dr52, (34a) 

"si? — JsRClsp dsR — JsRO'sp 

drVp = -fsRdrVR - fspdrVP " fssdrh, (34b) 
2 

VT 



I ~ivr2 \-'^sR^R + dsPVP + dspTip] , (34c) 
K[q — IJi 2 

8^62 = -esRtiR - espr]p - Csptip, (34d) 

at 0" = 1. The coefficients are given in (IC.llb -c) of Appendix O 

Two additional boundary conditions are imposed at the contact discontinuity. The requirement that 
there be no flow across the contact surface, that is u — dvc/dt = 0, implies drV^Src + 5vr — d^5rc = on 
each side of that surface. Upon substitution of the unperturbed solution we obtain 

dr5c = ^^i^^ _ (^ + 1) [1 + ^^(9^ In 5),] 5c, (35a) 
dr^c = qriR{T, cTc) - (m + 1) [1 + ac{da In G)c] 5c, (35b) 

where subscript c refers to values at the contact. Pressure balance across the contact discontinuity yields 

^p(r, Oc) - ip{T, Xc) = 2(m + 1) [a^d^ In F), - Xc{d^ In f)c]5c. (35c) 

After some manipulation of (|35b -c) we finally arrive at 

«9tCr(t, Xc) - 2qdrVR{T, (Jc) = -(m + 1)[1 + a^da lnG)c]£.R{T, Xc) 

+2g(m+ 1)[1 +Xc(9xln5')c]??/j(r,crc), (36a) 
r]p{T,ac) - ^p{t,Xc) =2{m + l){acdalnF)c- Xcdx'^^f)c}Sc, (36b) 



with 



. . N ^RjT, Xc) - 2qr]R{T, a^) , ^ 

'^^^ 2(m + l)[xc(5xln<7)c-^c(5.1nG)e]" ^ 



Equations (p3k -c). (p4k -d) and ([36k -c) provide a set of eight boundary conditions for the perturbation 
equations. 
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4.3 Numerical scheme 

To integrate equations ()25l) and ()27l) we first transform to a new set of variables, the so called "Riemann 
invariants" . In region 1 those are related to the old variables through 

e± = -^Cr ± ^Cp, (37a) 



6 = - - Hp/4., 



(37b,c) 



and in region 2 through 



V± 



m ^ 



1 



2qv«;7 



(38a) 



m = r]p - im-vp/i: m = m: (38b,c) 

with k(o") = W20/P20720' defined in Appendix IA.2I The equations for the new variables can be obtained 
upon appropriate transformations of ()25p . (I27p . and are derived in Appendix |Dj As shown there, 
propagate from the forward shock to the contact discontinuity while propagates in the opposite direction. 
Likewise, r]-^,r]3,rj4 propagate from the reverse shock to the contact discontinuity whereas 77- propagates 
from the contact discontinuity to the reverse shock. Thus, the two boundary conditions at the contact 
discontinuity, (136b . b). are applied to ^+ and r/_ after appropriate transformation of variables, and the 
boundary conditions at the forward and reverse shocks to the remaining Riemann invariants. 

Numerical simulations of equations (|D.3p and (ID.5|) were performed using the forward Euler scheme with 
upwind differencing of spatial derivatives. We also used for comparison a 4th order Runge-Kutta routine 
with an adaptive step size to advance the equations in time and found unnoticeable differences between the 
two methods for all cases studied below. For each experiment we have made several runs with increasing 
grid resolution until the result converged. The divergence of the unperturbed density at the contact have 
caused no difficulties. To test the code we exploited the analytic solutions for the spherical modes obtained 
below. We generally found excellent agreement. As an example, for the choice n = 1.1, k = 2 for which 
m = 0.645, and with 10^ grid points on each side of the contact discontinuity all of the Riemann invariants 
followed the analytic solution to an accuracy better than 10~^ up to a time r = 5 (or t/tQ = e^). The 
accuracy of the initial condition, specifically the accuracy at which the boundary conditions at the contact 
at r = were matched was 4.8 x 10"^ in this run. 

As a second test we solved ()D.3P for an impulsive BMK solution with k = 0. For this solution a contact 
surface does not exist {gx = 1 for every x)i so the boundary condition for ^_|_ at the contact needs to be 
replaced. We verified that the solution depends weakly on this condition provided it is fixed suffici ently 
far from the forward shock, at x ^ 1- We compared our result with that obtained in iGruzinovl (|200d ) and 
found excellent agreement!^. 



4.4 Spherical perturbations 

For spherical perturbations (/ = 0) the tangential velocity vanishes, viz., 5vt = 0, as can be directly 
seen from (|24li) and (f26l i). Then (f25]l . ([27l) admit solutions of the form (,aiXjT) = Caix)^'^^ j Vai^jT) = 
r]a{a)e'^'^; a = {R,P,p), where s is an eigenvalue determined from the boundary conditions at the contact 
discontinuity. There are two spherical modes. One is associated with a linear time translation of the 
unperturbed solution; that is, it is the difference between a solution at time t and a solution at time 



^ There is a small error i n Gruzinov ||200C ) . The coefficient of the 3rd term on the right-hand side of his 30 should be -17/3 (see l lB.3t )) 
and not -14/3 as in the astro-ph/0012364 version. However, we have confirmed that the effect of this error on the result is insignificant. 



12, 2009 7:53 Geophysical and Astrophysical Fluid Dynamics Levinson edit 
12 

t + 5t. Since to lowest order all surfaces propagate at the speed of light the distance between corresponding 
surfaces of the two solutions should remain 6t at all times. This implies 6Va = to order 0(r~^) and 
from ()29p we anticipate that s = —{m + 1) for this mode. From dH) and (I16p it can be readily seen 
that the change in the self-similarity parameters x and a corresponding to a time translation 5t is 5x = 
2(m + l)T\5t/t = 2(m + l)Be-(™+i)^(5t and Sa = 8x1x2- To order 0(rj-2) equations ((EJi-c)) yield 
= {dxli)8x-, 8pi = {dxP^o)^X^ and 6p'i = {d^p'iQ)5x- Note that since ^xd^ = —Stdr (see (jA.lb )) 
the latter relations simply mean that the Lagrange perturbations of the flow parameters vanish; that is, 
Api = 5pi + dtdrPio = Spi — {d^pio)6x = 0, etc. By employing ([2lb -c) we finally obtain 

^r{x) -- 
Ux) -- 
Ux) -- 

in region 1. Likewise, in region 2 we find, using 

mix) = 

Vpix) = 
Vpix) = 

By direct substitution it can be shown that 
Caix) and Tjaix) given by ([39]) and (jlO]) . respectively, satisfy the perturbation equations ([25|) . ([27|) . and all 
the boundary conditions. Note that for this mode 6r2Q = 6r2 in ()C.3b ) and ()C.7b ). 

The second eigenmode of order / = was found by numerically integrating the set of ODEs obtained 
upon substitution of the relations ^aiXi^) = ^a{x)^^^ ^ i]a{o',T) = r]a{cr)e^'^ into equations (j25]) and (j27|) . 
The method of solution was to guess the eigenvalue s and integrate the equations in each region from 
the shock to the contact surface. The process was repeated until the condition ([36b ) was satisfied to the 
required accuracy. For each choice of parameters we found, using this method, two solutions; the first one 
coincides with the analytic solution given by ([39|) and (|10|) . The second one has a smaller decay rate, 
s > — (m + 1). 

4.5 Non- spherical perturbations 

For the simulations of non-spherical perturbations (/ 7^ 0) we employed the analytic solution given in (j39p - 
()40p as the initial condition. We have also made some runs for comparison using the second spherical mode 
found above as the initial condition and verified that the long term evolution of non-spherical growing 
modes is essentially independent of the details of the initial state. The free parameters of the model are: 
the indices k and n characterizing the density profiles of the unshocked ambient medium and ejecta (see 
figured]), and the adiabatic index 7 of the shocked ejecta. The exponent m and the ratio q are given by 
(|2Up and (|15p . respectively, for a given choice of n, k and 7. Since the reverse shock is non or at best mildly 
relativistic in the cases examined below we adopt 7 = 5/3. We have explored solutions for a range of values 
of n and k for which m lies in the range 0.5 — 2. The results displayed in figures [3]|6] were computed using 
our canonical choice of parameters: n = 1.1, k = 2, for which m = 0.645, q = 1.06. The corresponding 
unperturbed solution is displayed in figure [2] (left panel). For other values of n and k the solution exhibits 
the same qualitative behavior, albeit with a larger growth rate for larger values of m . The dependence of 
the convective growth rate on m is examined in figure [7] and discussed below. 

Quite generally we find that eigenmodes having angular scales larger than the causality scale, roughly 
l/Ti < 1, decay with time. This is expected since for these modes information can only propagate over 
distances smaller than the characteristic wavelength, so that a causal section evolves like a spherical mode. 
An example is shown in figure O where solutions obtained for 1(1 + l)/r^Q = 0.2 (upper panels) and 



:2{m + l)d^lng, (39a) 

-.2{m + l)d^lnh, (39b) 

: 2(m + l)a;^ln/, (39c) 

dHJi-c) 

^"^^^^ d^lnG, (40a) 

2{m + l)d^\nH, (40b) 

2{m + l)d^ln F. (40c) 



{r,x) = ea(x)e-(™+^)^ r/„(r,x) = ^?a(x)e^("+')^ with 
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Figure 3. Time evolution of the perturbations for n = 1.1, k = 2, l{l + l) /Fjq =0.2 (upper panels) and + l)/T-fg = 1.5 (lower panels). 
Here Fio being the initial Lorentz factor of the forward shock. The dimensionless distortions of the different surfaces are delineated in the 
left panels and the relative pressure perturbations (i.e., normalized to their initial values) in the right panels, as indicated. The dashed 
red line corresponds to the analytic solution of the spherical mode and is plotted here for a comparison. 

+ l)/r^Q = 1.5 (bottom panels) are exhibited. Here Fio = ri(r = 0) denotes the initial Lorentz factor 
of the forward shock. The analytic solution for the spherical mode is plotted for a comparison in the 
right panels (dashed red line). The evolution of the dimensionless distortions of the forward and reverse 
shock fronts, (5?"i(2)/'''i(2) = ^U2)/'^i(2)^ contact surface, brdvc = bd^^^ (see equation (p8]) ) are 

displayed in the left panels. Relative pressure perturbations are shown in the right panels. The decay of 
perturbations with an angular scale sufficiently larger than the initial causality scale is evident from this 
example. The deviation from the analytic solution is small, as seen in the upper right panel. Evolution of 
longer wavelength modes, lij, + l)/r^Q < 0.2, is practically identical to that of the spherical mode. The 
turnover at r ~ 1.5 exhibited by the solution displayed in the bottom panels is due the fact that /(/ + l)/r^ 
increases with time as exp(mr), so that the wave enters the 'horizon' in this case early enough to affect the 
evolution of the perturbations at the contact surface. This is roughly the boarder case separating stable 
and unstable modes. 

Modes of order / > Fio are found to be unstable. This is demonstrated in figure [H where the evolution of 
an eigenmode of order I = lO^Fio is delineated. Oscillations resulting from sound waves crossing, followed 
by a rapid growth of perturbations in the shocked ejecta and near the contact in the shocked ambient 
medium are clearly seen. The frequency of oscillations is found to be proportional to the spherical harmonic 
degree I, as anticipated. Moreover, the frequency of oscillations in region 2 is smaller than that in region 1, 
owing to the difference in specific enthalpies of the fluids on each side of the contact surface. The growth 
of perturbations at the contact surface commences very early on, as seen in figure HI The instability then 
propagates from the contact to the forward and reverse shocks. The reverse shock responds rather quickly, 
mainly because it is located much closer to the contact than the forward shock. The propagation of the 
instability in the shocked ambient medium (region 1) is rather slow, as illustrated in figure [5l The reason 
is that energy pumped from the contact into this region is now distributed over a much larger volume. 
However, in reality the instability near the contact will quickly reach the nonlinear regime and saturate, 
at which point the linear analysis breaks down. What then anticipated is formation of R-T fingers that 
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Figure 4. Same as figure 2 but for 1(1 + l)/r^Q = 10*. Left panel: dimensionless distortions of the forward (dotted-dashed line) and 
reverse (dashed line) shock fronts and the contact surface (solid line). Right panel: relative pressure perturbations at the forward shock 
front (dotted-dashed line), reverse shock front (dashed line), and the contact surface (solid line). 



expa nd to the forward shock, as seen in the nonrelativistic case (jChevaUer et a/.lll992l . lJun and Norman 
19961 ). To study the response of the forward shock requires fuh, high resolution MHD simulations. 



The growth of the perturbations that starts following the decay of the transient initial state is exponential 
in the time r, or algebraic in the physical time t, viz. 5Q oc e**"^ oc {t/toY where Q represents any of the fluid 
quantities. This is shown in figure [6] for the relativistic distortion 6c defined in (j28p . Growth of variables 
located at larger distances from the contact commences at later times, but follow with the roughly same 
growth rate. It is also evident from figure [6] that the growth rate s increase with increasing mode degree 
I. From our numerical simulations we find the scaling s oc y^Z/Fio, with the proportionality constant 
depending predominantly on m. Examples are exhibited in figure [7] for two different cases; the first one 
corresponds to m = 0.645 (solid line) and the second one to m = 2 (dashed line). A fit to the lines in fig [7] 
gives s ~ y/l/^io in the former case and s ~ 2.45y //Fio in the latter. This scaling is expected for a R-T 
instability which is driven by the "effective" gravitational force felt by the decelerating contact interface. 
To see this, let Oc denotes the 4-acceleration of the contact interface and Uc its 4- velocity. In flat spacetime 
the acceleration is orthogonal to the velocity, a^Uf^ = 0. In the rest frame of the contact = (0, a^,0, 0) 
from which we obtain = (a^)^. From these relations and the normalization of the 4- velocity, = 1, 
we get a'j. = du^/dt ~ dFc/dt = -mTc/2t for F^ oc (t/to)"™"- The "effective" gravitational force felt by 
the decelerating contact is just g = —a!^. Now, the wavevector component parallel to the contact surface, 
/c||, is given roughly by = rdl for an eigenmode of order /. We thus have gk\\t' ~ \J ml/2T(. oc s, 
noting that t' = t/Tc and that the ratio Fc/Fi is independent of time. Consequently, the growth rate, as 
measured in the rest frame of the contact, is s/t' oc y/gk\\, which is just the R-T growth rate in the case 
P2c S> pic- Note that the dependence of s on m, as derived from the simulations, reflects not only the effect 
of deceleration but also the dependence of the unperturbed solution on m, in particular the density ratio 
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Figure 5. Profile of the dimensionless perturbation of radial velocity at different times. Left panel: shocked ambient medium (region 1). 
The location of the forward shock is at x = 1- The location of the contact discontinuity Xc is indicated. Rapid growth of the perturbation 
near the contact is clearly seen. The wave propagating from the forward shock to the contact interface corresponds to transmission of 
the initial shock disturbance. Right panel: shocked ejecta (region 2). Here riji is plotted against the normalized variable cr = x/X2 (see 
equation jlOjl). The locations of the reverse shock and the contact discontinuity are at tr = 1 and a = cTc = Xc/X2, respectively. As seen, 
at T = 1 the instability already propagated throughout the entire region. 
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Figure 6. Absolute value of the relativistic distortion 5c defined in equation l|28p (see also l(36j;)) versus time t, for different values of 
the spherical harmonic degree /. 



at the contact, as can be seen from figure [2j 
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Figure 7. Dimensionless growth rate versus (/Fiq. The solid line corresponds to m = 0.645 and the dashed line to m = 2. 

5 Implications for GRBs 

In the standard GRB model the afterglow emission observed following the ex plosion is commonly attri buted 
to s ynchrotron cooling of relativis t ic electrons behiii d the forward shock (M eszaros and Rees 19971 ) (but 
c.f., ( Lucas and Beloborodov 2007 . Genet et al. 20071 )). The observations seem to indicate the presence of 
strong magnetic fields over a significant portion of the shocked circumburst layer. The magnetic energy 
density estimated from the d ata, roughly a fraction ~ 10~^ — 10~^ of the internal energy density 
( Panaitescu and Kumar 20021 ). is several orders of magnitudes larger than that expected from compression 
of the preshock magnetic field. This implies magnetic field generation or amplification by some mechanism. 
Whether kinetic processes can generate such magnetic fields is yet an open issue. Plasma instabilities that 
develop in the collisionle ss shock transitio n generate strong magnetic fields on kinetic scales. However, 
recent shock simulations ( Spitkovskv 20081 ) indicate that the fields thereby produced decay rapidly over a 
few skin depths, before reaching the MHD scale. Magnetic field generation by streaming ultra-relativistic 
protons (or cosmic rays) in the immediate shock upstream has also been considered. In this scenario the 
magnetized cosmic ray precursor is envisioned to form as an inherent part of the shock transition. The 
key question here is whether the shock can indeed inject enough energy in the form of ultra-relativistic 
particles to sustain the required large scale magnetic field. To study this process using self-consistent shock 
simulations requires computing time beyond present capabilities. 

An alte rnative to plasma instabilit ies is magnetic field amplification by MHD turbulence . It has been 
proposed (jSironi and Goodmanl[2007l ) that macroscopic turbulence might be produced via the interaction 
of the forward shock with a clumpy circumburst medium. In this scenario the response of the shock to 
the preshock de nsity inhomqgenei ties leads to generation of vorticity and the consequent amplification of 
magnetic fields ( Zhang et al. 20091 ) . Amplification of magnetic energy to the level inferred from observations 
requires large (order unity) density contrasts and filling factors of the clumps. Whether such conditions 
exist in the surrounding environment is yet an open issue. Here we propose that the convective instability 
found above may be an inherent source of turbulence in the shocked circumburst layer, at least at early 
times. The instability may also lead to nonlinear distortions of the shock front itself without the need for an 
external driver. If the ejecta is magnetized at a level smaller than that required to suppress the instability 
but still much larger than that of the unshocked ambient medium, then mixing of the magnetized ejecta 
with the shocked ambient gas via growth of R-T fingers alone can lead to strong magnetization of shocked 
cirumburst layer at sufficiently early times. How the system evolves at later times, after the reverse shock 
crosse s the ejecta, is unclear at present. The stability analysis of the BMK solution performed by Gruzinov 
suggests that it may be a very slow attractor. Linear perturbations of the forward shock in the 
BMK phase decay very slowly. Whether this behavior persists also in the nonlinear regime is yet to be 
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determined. If it does then it could well be that the growth of R-T fingers and, perhaps, the nonlinear 
oscillations of the forward shock itself which are induced by the co nvective instability may be a source 
of vorticity during a long portion of t he evo lution of the blast wave (iMilosavljevic and Nakail [2007! ^ ■ The 



simulations performed by Zhang et al. ( 20091 ) suggest a saturation level of ~ 5 x 10 ^ for the turbulence 



induced magnetic energy density, weakly dependent on the initial magnetic field strength. 

The convective instability may also affect the emission from the shocked ejecta. As shown above, the 
effect of the instability on the reverse shock is prompt and dramatic. The nonlinear distortions of the 
reverse shock may strongly alter particle acceleration and the emission processes during reverse shock 
crossing. It is tempting to speculate that the lack of observed optical flashes, that are anticipated in the 
"standard" model, and the behavior of the early afterglow phase may be attributed to the instability 
discussed here, although we do not offer at present any specific explanation. At any rate, the salient 
lesson is that a careful analysis that takes account of this process is required to better understand the 
observational characteristics of the emission during the early post-prompt phase. 



6 Summary 

We have performed a global linear stability analysis of a self-similar solution describing the interaction of 
a relativistic shell with an ambient medium. Our analysis indicates a strong convective instability at early 
stages of the evolution of the dense ejecta as it sweeps a lighter ambient gas. Our main findings are: 

1. Eigenmodes having angular scales smaller than the causality scale, roughly l/Ti > 1, where Fi is the 
Lorentz factor of the blast wave, are unstable and exhibit a rapid growth. Lower order modes for which 
l/Ti < 1 are stable. 

2. Growth of perturbations starts promptly near the contact discontinuity. The instability then prop- 
agates towards the forward and reverse shocks. The reverse shock responds quickly to the growth of 
distortions at the contact. Propagation of the signal to the forward shock is much slower. The instabil- 
ity near the contact becomes nonlinear well before the signal arrives at the forward shock, so full MHD 
simulations are needed to resolve the effect of the instability on the forward shock. 

3. The growth is algebraic in time, that is, SQ oc {t/toY fluid quantity Q. The dimensionless 
growth rate scales as s oc yU/Ti, with a proportionality constant that increases with increasing m. This 
implies development of a very small scale structure with a significant amplitude, up to the dissipation scale, 
at least at sufficiently early times. The effect of such corrugations on the collisionless shock transition and 
related processes, particularly particle acceleration, needs to be explored. 

Unfortunately, the linear analysis outlined above is restricted to a li mited set of conditions under which 
the unperturbed self-similar solution of Nakamura &; Shigeyama ( 20061 ) is applicable. It is naively expected 



that the instability will be strongly suppressed in cases where the ejecta is highly magnetized and/or if 
the reverse shock is highly relativistic. Full 3D MHD simulations should be exploited to study this process 
in other situations, and to follow the evolution of the convective instability in the nonlinear regime. As 
illustrated above, high resolution simulations that can resolve angular scales A9 « 1/F are required, 
posing a great numerical challenge. We believe that our findings strongly motivate such efforts. 

I thank A. Ditkowski for a technical help in the development of the code, and M. Aloy, D. Kushnir, M. 
Lyutikov, A. MacFadyen, E. Nakar and E. Waxman for enlightening discussions. This work was supported 
by an ISF grant for the Israeli Center for High Energy Astrophysics, and by the NORDITA program on 
Physics of relativistic fiows. 
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Appendix A 

Derivation of the unperturbed flow equations 
A.l Region 1 



The derivation of the flow equations i n regio n 1 (the shocked ambient medium) is straightforward and 
follows that in Blandford and McKed ( 19761 ). Transforming from the coordinates {r,t) to {Xit)^ with 



T = lnt and x given by ^ , and using the relations 

tdt = dr + [{m + l)(2r2 -x) + l]d^, (A.la) 

tdr = -[I + 2{m + l)Tl]d^, (A.lb) 

id° = dr + {m + l){2/g - x)d^, (A.lc) 

one obtains, upon substitution of equations (fT2b -c) into the flow equations (l3^,b), 

1 din (7 (7m + 3/c — 4) — (m + 2)gx 



g dx {m + l){4-8gx + g^X^) 

1 d In /i _ 2(9m + 5A; - 8) - 2(5m + Ak - 6)gx + {m + k - 2)g'^x^ 
g dx {m + l){A-Sgx + g'^X^){2- gx) 

1 d In / _ 8(m - I) + Ak - {m + k - 4)gx 



g dx {m + l){4-8gx + g'^X 

A. 2 Region 2 

In region 2 we use the coordinates (r, a). Then, 



(A. 2a) 
(A.2b) 
(A.2c) 



tdt = dr + [{m + l){2Tl -a) + l]d^, (A.3a) 
tdr = -[l + 2{m + l)Tl]da, (A.3b) 
td'l = dr + {m + l){l/qG-a)d^. (A.3c) 



Substituting p8b -c) into equations (l3^,b) yields 

^ ,dlnF ^ , dlnG (mn — Km — 6) ^ /. . n 

2{l + qGa)— {l-qGa)K—— = ^ -— '-qG, A.4a 

do- dfj [m + 1) 

2 1-(zG(7— 2-— = - ^ ' qG, A.4b 

da da [m + 1) 

dlnF dlnG [7(m + 2) - mn - 6(7 - 1)] /a . n 

2 1-gGa— j{l+qGa)— = ^ ' '-^qG, A.4c 

da da [m + 1) 

where K{a) = W^jvil^^ and 7 denotes the adiabatic index. 
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Appendix B 

Derivation of the equations for the dimensionless perturbations 
B.l Region 1 

To derive the equations for the dimensionless perturbations {t, x) first write equations (j23b -d) in the 
new coordinates (X)''")) using the relations ()A.lb -c) and (|24b-dl. Recalling that Wi = ^pi^i and noting 
that rV • SviT = + Ij^T^im; we obtain from the continuity equation (j23h ) 

2(m -I- 1 2(m + l^l 
dAp + {m + l){2/g - x)d^ip - ^ ^ ' + ^ ^ Ing-d^ In h)(n = 1(1 + (B.la) 



Equation (p3b ) gives 



2 (m -\- ^'] 1 



+^-^[-9^1n5 - d^\^i\iR = -l{l + IKt, (B.lb) 

the transverse component of the momentum equation ()23b )l gives 

dr^T + {m + l){2/g - x)d^CT + [{m + l){2/g - x)d^ Hfg/h) - m]^T = (B.lc) 

and the energy equation ([23l i) gives 

dri^p + 2Cr) + 2{m + l){2/g - x)d^CR - {m + l){2/g + x)d^Cp + CoR^r + CoP^p = 0, (B.ld) 

where 

CoR = -{3m + k)-{m+ l)xd^[2{\ng) + (In /)], (B.2a) 

CoP = -(3m + k)+ 2(m + l){2/g - x)5^(ln <7) -{m+ l)(2/5 + x)5x(ln/). (B.2b) 

Note that (|RTb l holds only for 0. In this case it is readily seen from (jB.lb .c) that S^p and must 
have different time evolution, owing to the extra factor on the right-hand side of (jB.lb ). implying 
breakdown of self-similarity. In the spherical case (/ = 0), Vrl^^ = 0, and equation (j23b l is identically 
zero. The right-hand sides of (jB.lb .b) then vanishes allowing separation of variables. 

For an impulsive BMK76 solution with k = 0, m = 3, g = x~^} ^ = ^-7/4 j _ ^^-17/12 above 
equations reduce to 

drCp + ^xd^Cp - Sxd^^R + GCr = l{l + 1)Ct, (B.3a) 

dr[(p + (2/3)^] + 4x9x?P - 8x9;,^^ + 6U = ^1(1 + mr, (B.3b) 



3 2gT^ 



drir + ^xdxiT - -er + 77^ = 0, (B.3c) 



14 

driip + 2^r) + Sxd^^R - 12xd^^p + = 0, (B.3d) 



as derived originally by Gruzinov ( 2000l ) (except for the 3rd term on the right-hand side of ()B.3b ). see 



footnote at the end of section [^3]) . After some manipulation of equations ()B.lh -d) we arrive at (j25p . with 



12, 2009 7:53 Geophysical and Astrophysical Fluid Dynamics Levinson edit 
20 

the different coefficients given by 



and 



Arr = App = {m + l){x-^/g), (B.4a) 
Apr = \Arp = 2ApR = l(!!i±l), (B.4b) 
App = ATT = -{m + 1) (2/5 - x) , (B.4c) 



3(3m + A;) (m + l)/2 3(m + l)/2 , , 

Brr = ^ + ^^-^ f - + 3xJ ln5 - ^ ^ ^ U " ^ ' ^^-^^^ 

Bpn = -^^^ - (m + 1) + x) 5xln. + ^ - X ) 9, In/, (B.5b) 



3 i^m + k) , ^ /2 \ „ (m + 1) /2 , , , 

Bpp = --Brp = ^ ^ + (m + 1) - xj 5^1n<7 - ^^-^ (^- + xj In/, (B.5c) 

= -2S^r = 2BpT = 2l{l + 1), (B.5d) 

BpR = -^^d^Hg/h), (B.5e) 

Btp = -l/{2gr^), (B.5f) 

= m - (m + 1) - 9x Hfd/h). (B.5g) 



All other coefficients vanish. 



B.2 Region 2 

The derivation of the perturbation equations in region 2 is similar to that in region 1, but slightly more 
involved. Here we express (f23k -d) in the coordinates (o", r), using the relations ()A.3b -c) and (f26k -d). To 
shorten the notation we define a = 7/(7 — 1) and A± = (m + l)(l/gG it a). From the continuity equation 
(I23h ) we then obtain 

2(m -\- 1 2(in -\- 1 
drVp + A^d^7]p - ^ J ' d^m + ^ ^ \ d„ InC - 5, Ini/)??^ = /(/ + (B.6a) 



noting that rV • (5v2t = + l)fyT^im- Likewise, ([25b ) yields 



9^r/T + A-daTlT + 



-A_d„lniF/H) + ^^i^A_9<,lnG-m/2 
K 2k 



where the relation W20/P20 = KO'CFg has been used, the transverse component of the momentum equation 
(I23b) gives 

drivP + iqm) + ^-daVP - iq^+daVR + Q ' [^d,\nG- 2a^lnF]7?R = 7/(/ + l)r?T, (B.6c) 
and (f23]i) gives 

nqdrfiR + dr-qp + KqA^d^ilR - A+d^tlP + Qi^r/R + Cipijp + QpT/p = 0, (B.6d) 
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with 



K — a 



CiR 



2(m + 1 



G 



(ui _|_ X) 

dalnF - K — — — d„ In G - qcip. 



(B.7a) 
(B.7b) 



After some algebra equations ()B.6b -d) can be recast into the form of (j27p . with the coefficients Cap, 
Daf} given exphcitly by: 



and 



C 



RR — Cpp 



{m + 1) 



K + 7 



Crp 



CpR 



qG{K-^) 
2(m+ 1) 



C, 



pR 



2(m + 1) 
G 



Drr 

DpR 

Dpp 
DpT 
DpR 



Cpp = CTT = {m+l){a-l/qG), 



""^^^ [-Ad^lnF + {k + j)d^lnG] + — ^ [A_5^ hi G-m] 



qG{K-^) 
2(m + 1 



2(k-7)' 



G(^-7) 
-qjDpp 



[(7 + K)d^ In F - K75^ In G] - ^^^"^ [A_9^ In G - m] , 



2(k-7) 



-qKDpT 

2(m + 1) 
G 

£'pT = + 1), 

Dtp 
Dtt 



- a) 
2(^^-7) 

7^^(/ + 1) 
(k -7)g ' 

5,ln(G/F), 



A_a^lnG-m], 



KqVlG' 

--A_daln{F/H) 



2k 



;A_a^lnG + m/2. 



(B.8a) 

(B.8b) 

(B.8c) 
(B.8d) 

(B.9a) 

(B.9b) 

(B.9c) 

(B.9d) 

(B.9e) 
(B.9f) 
(B.9g) 

(B.9h) 



Appendix C 

Boundary conditions at the reverse shock 

The normal to the reverse shock is written = + 6n2p, where 71-2^ and Sn2^ are given by equations 
(fSTk .b). respectively, with the subscript 1 replaced by 2. Applying (f30k ) to the reverse shock we have 

P>20^n2p + {A2p'2vl^o + p'2oA2<)nO = p',v^6n2p + (Aap'^^ + /o^AaOnO, (C.l) 

with the density of the unshocked ejecta pe is given by ([7]), Ve = (1, fe, 0, 0), etc. The operator A2 is defined 
below equation (I29p . Now, the velocity of the unshocked ejecta near the shock surface evolves with time 
as dve/dt = dtVe + (V2 + 5V2)drVe- Integrating over time from t = to we find A2Ve = {6r2/t){l — 5r2o/Sr2) 
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to the order to which we are working, where 6r2o is the initial displacement of the reverse shock surface. 
From (f28l) we then obtain A2fe?^2/i = (1 ~ '5?^2o/'^?'2)(<^2/r2)5^/m- The corresponding convective change 
of the density p'^ is A2lnpg = (1 — n)(m + l)r2A2fe- Using ([3Tb) with the subscript 1 replaced by 
2 we also have 2(m + I)ve5n2^ = —{m + 2)T25V2Yim. For the perturbations of the shocked ejecta we 
obtain, using the unperturbed flow parameters (fT8h -c). A2ln/92 = [vp — 2(m + l)523(j(lnff)2]llr«; and 
gr2A2t'2 = [qiIR — im + l)(52i9cr(ln G)2]l/m- Substituting the above results into (jC.ip yields 



where 



drS2 + AsRm + AspTJp + Ass62 = 0, (C.2) 



AsR = -^Asp = ■ (C.3a) 

q — I q — m — I 

Ass = m + 1- . ^ ^ j a.(lnG)2 + - l)5.(lnF)2} 
(g — m — 1) 

(m + l)(m + 2 - mn)(g - 1) rr^ QK^ 

^7 7^ [I- 6r2o/6r2). (C.3b) 

2[q — m — 1) 



The transverse component of (j30b ) reduces to 

= - 1) Fi 

and the zeroth component reads 



W20vt^QSn2p + (A2VF2'V20 + W^20A2V2 )"2/i " (P20'^f^20 + ^2P2n2o) = 

WeV^Sn2p + {A2WeV^ + ^eA2^;^)n2^, (C.5) 

with We = Pell A2We = Pell{2 - n){m + 1)52(1 " 5r2o / 5r2)Yi^, A2W2 = P2o{k + a)q'^iTlA2V2) + qin - 
a)A2 lnp2 + qcLA2 lnp2, and A2 lnp2 = [^p — 2(m + l)62da{lia -F)2]^m- Rearranging terms we finally arrive 
at 

BstdrS2 + BssS2 + BsRVR + Bgprjp + B^piip = 0, (C.6) 



where 



Bst = l ;^^+ o D ' (C-7a) 



S.i? = ^^^^(ac - a) + aq\ (C.7b) 
i^..^^^^^, (C.7C) 
i?.p = ^('^-a), (C.7d) 



Bs5 = (m + 1) 



- :^a,(lnG) - 2Bspd^{\nF) - 2Bspd^{\nH) 



+-^{m + 1) (mn/2 - m - 1) (1 - 6r2o/6r2). (C.7e) 

-^20 
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Subtracting the radial component of (130b ) from (ICSh we obtain 

[W2ov^o5n2^, + (AaVFa^o + W2oA20n^^](t;2o - Ve) + t^2o(A2^;2)^^^o"2M 



+P20(V2 - Ve)rl6V2 + A2P2(1 " ^2t'e)r2 = We{A2VeKnL, (Ci 



whereby the final condition, 



Cstdr52 + CsSh + CsRVR + CsPVP + CspTJp = 0, 



is derived, with 



Cst 

CsR 
CsP 

Can 



(m + 1 - q){q + 1)k 
{q — m — 1] 



m, 



[{q + 1)k + iq- l)a] + (q - l)(m + 1)k, 



{q-m- l){q - 1) 
2q 

{q-m-l){q + l) 



2q 



a + (m + 2). 



{k - a), 



CsS = (m + 1) 



C. 



c. 



sR 



St 



d^ilnG) - 2Cspd^{lnF) - 2C,pd^{\nH) 



(C.9) 

(C.lOa) 
(C.lOb) 
(C.lOc) 
(C.lOd) 
(C.lOe) 



With the exception of the spherical mode associated with a linear time translation of the solution outlined 
in section [3l for which 5V2 = and 6r2{t) = 5r2Q (see section we shall assume 5r2Q = 0. We can then 
eliminate 82 from ()C.2p . ()C.4p . ()C.6P and ()C.9P to obtain equations ([Mb -d). where we define 



fs\ 



CsaBst ~ BsaCst _ 

CssBst — BgsCst ' 

CsaBsS BgQ^CgS ^ 

CssBst — BssCst 
Bsx — AsxBst ^ 
Bsp ~ AgpBgt 



a = {R,P,p), 



a = {R,P,p), 



X = {R,P,6). 



(C.lla) 



(C.llb) 



(C.llc) 



Appendix D 

Equations for the Riemann invariants 

Define the vector |^) = {(,r,(,p,S,pi(,t)- Then ([2^ can be expressed in matrix notation as 

dr\0=Ad^\0 + B\0, 



(D.l) 



with A{x) and B{x) denoting the matrices (Aa/s) and (Ba/s), respectively. Let {ipg); q = (+, — , 3, 4) be the 
eigenvectors of A and Xg the corresponding eigenvalues. We find 



A± = Arr ± ^ApR, 
A3 = A4 = -(m + l)(2/(7-x), 



(D.2a) 

(D.2b) 
(D.2c) 
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where ^_rr and Ap^ are given by ()B.4b .b) respectively. The Riemann invariants in region 1, given exphcitly 
in ([37ll . are the coefficients of the expansion of |^) in the basis vectors iV'g), that is, |^) = S^^gl^g). They 
obey the equations 

where Rqp = (^g|-B|^/^p). We find that A_, A3, A4 are negative everywhere in region 1, whereas A-|- is positive, 
so that ^-,^3,^4 propagate from the forward shock to the contact while ^+ propagates in the opposite 
direction. 

The analysis in region 2 is similar but slightly more complicated, as the eigenvectors \C,p) > of the matrix 
C, the elements of which are given in ()B.8b -d). depend on the self-similarity parameter a. The Riemann 
invariants defined in ([38|) were computed using the relation |r/) = T,qr]g\(g), where |r/) = {r]R,rip,rip,r]T)- 
For the eigenvalues of the matrix C we have 

~X± = Crr±-^^, (D.4a) 
A3 = A4 = -(m + l)(l/gG-a). (D.4b) 
The equation governing the evolution of the Riemann invariants in this region reads: 

drVq = ^qdxVq + ^pTqp'Hp, (D.5) 

where Tqp = {C,q\D\Cp) + {C,q\C\daC,p) ■ The eigenvalues A+, A3, A4 are positive everywhere while A_ is negative 
hence ?7+ , 773 , 774 propagate from the reverse shock to the contact while rj- propagates from the contact to 
the reverse shock. 
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